Analytical expressions for the Electromagnetic 
Dyadic Green's Function in Graphene and thin 

layers 

A. Yu. Nikitin, F. J. Garcia- Vidal, and L. Martin-Moreno 



Abstract — An analytical general analysis of the electromagnetic 
Dyadic Green's Function for two-dimensional sheet (or a very 
thin film) is presented, with an emphasis on on the case of 
graphene. A modified steepest descent treatment of the fields 
from a point dipole given in the form of Sommerfeld integrals 
is performed. We sequentially derive the expressions for both 
out-of-plane and in-plane fields of both polarizations. It is shown 
that the analytical approximation provided is very precise in 
a wide range of distances from a point source, down to a 
deep subwavelength region (1/100 of wavelength). We separate 
the contribution from the pole, the branch point and discuss 
their interference. The asymptotic expressions for the fields are 
composed of the plasmon, Norton wave and the components 
corresponding to free space. 

Index Terms — Graphene, thin films, plasmon, Dyadic Green's 
Function. 

I. Introduction 

ELECTROMAGNETIC properties of graphene have re- 
cently received a lot of attention due to a variety of 
application in photonics (TJ. One of the attractive properties 
of graphene is its capacity to support highly-localized (nano- 
metric) surface modes, i.e. graphene surface plasmons (GSPs) 
in terahertz (THz) and micro- wave frequency ranges p|-p2|. 
GSPs and their gate tunability open interesting possibilities 
for construction of tunable meta-materials (4), and merging 
photonics and electronics. In particular, for applications re- 
lated to strong light-matter interactions and biosensing, the 
interaction of a graphene sheet with a point emitter presents 
a special interest (5|-|T2|. Localized excitation has allowed 
experimental demonstration of GSPs (TTJ, (12). In comparison 
with other experimental techniques, local excitation of GSPs 
is more favorable due to very high values of GSP momentums. 

The computation of patterns of the electromagnetic fields 
in graphene created by point source (as well as spontaneous 
emission rates), requires the knowledge of the Dyadic Green's 
Function (DGF) (3), (3), (5)-(T0). Its calculation involves 
notoriously difficult Sommerfeld- type integrals, with the in- 
tegrands containing quickly oscillating functions, poles and 
branch cuts p3l-p5|. 

In this paper we will perform an analytical analysis of DGF, 
providing an asymptotic series expansion with the modified 
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steepest-decent method (13), (15). We will show that this 
expansion, while being exact for long distances, in practice 
is very precise outside of its formal validity range. We will 
explicitly provide contributions for GSP, out-of-plane prop- 
agating waves and field components decaying algebraically 
along the graphene sheet. Notice that while all the examples 
correspond to the conductivity of graphene, the analytical 
expressions are applicable to any two-dimensional (2D) sheet. 
The analytical expression can be useful for treatment of 
more complicated problems, related for instance to Lippmann- 
Sch winger integral equation (14). 




Fig. 1. The schematic of the studied system. A point source (dipole) is placed 
onto a graphene sheet, characterized by a two-dimensional conductivity a. 
The colorplot represents an example of the spatial distribution of the electric 
field modulus for shown orientation of the dipole moment. Representative 
parameters for the graphene sheet corresponding to the formation of GSP are 
(see Appendix B for the definitions): fl = 0.4, T = 300 K, r = 1 ps, 
fj, = 0.2 eV. 



II. Formulation of the problem 

Let us consider a dipole with an arbitrary oriented dipole 
moment p placed at the point (0,0, Z'\ with Z' being the 
distance from a free-standing graphene sheet which covers the 
plane Z = 0, see Fig [I] Without loss of generality, we will 
suppose that Z' < 0. The time dependency is supposed to 
be e~ lUJt , where u is the angular frequency. Throughout the 
article, we express coordinates in the in-plane (R) and normal 
(Z) directions to the graphene sheet in dimensionless units as 
r = A'^.R and z = k u Z, with k u = 2tt/X being the free- space 
wavevector. 
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Graphene is represented by its in-plane complex conductiv- 
ity a. Throughout this paper, all examples are performed for 
a based on the random-phase-approximation p6|-|T8), see 
Appendix B. 

The electric field E(r, z) emitted by our electric dipole, is 
given through the DGF G(r, z, z') = G(r, z\ r' = 0, z') by 
the following relation 



E(r^) = G(r,^')p(A 



(1) 



In the next section we will provide the exact general expres- 
sions for G(r, z; z'\ 

III. The general form of the Green's dyadic 

We have performed the analysis for a graphene placed onto 
the boundary of two different dielectrics. However, we have 
found that there is no much qualitative difference between this 
general case and free-standing (suspended) graphene. Since 
the analytical formula in the general case are quite lengthy, in 
this paper we will consider the DGF for a suspended graphene 
(taking also into account that free-standing samples are widely 
used in the experiments). 

DGF satisfy the following differential equation 

V x V x G(r, z\ z') - *£G(r, z\ z') = 16(r)6(z - z'), (2) 

where 1 is a diagonal unit matrix and S is the Dirac delta 
function. The DGF must be complemented by the boundary 
conditions at the graphene sheet that we present below. 

A. Angular representation of DGF in Cartesian coordinates 

The solution for G can be can be expressed (see e.g. (l3) , 
fl4] l) in terms of plane waves in free-space u clr e' lcir+tqzZ , 
characterized by their in-plane momentum q = k/fc^ and 
polarization r = TE, TM). The unitary vectors characterizing 
the polarization of each mode are: 

/ Qx 



A qTE 




U 



qTM 



q 



(3) 



where "+" ("— ") upward (downward) propagation along z 
respectively. q z = y/l — q 2 is the normalized z-component of 
the wavevectors. 
The DGF reads 

G(r, Z ] z f ) = G (r, z\ z f ) + G R (r, z\ z f ), z f < 0, z < 0, 

G(r, z\ z) = G T (r, z\ z), z' < 0, z > 0, 

(4) 



with Go being the DGF in free space (FS), 

G (r,,;,') = £/ g u ±„±r.^.l«-»'l 



(5) 



and Gr, Gt being the contributions due to the reflection and 
transmission in our 2D system 



G R (V,Z;Z) = ^2 j ^^ U qr U , 

G T (r,,;,0 = E/^ T ;<u, 



+T iqr-iq z (z+z') 
qr c i 



+T icp+iq z (z-z') 
qr e 



(6) 



In these expressions the superscript "T" stands for transposi- 
tion. 

In the above expressions R T q and TJ" are the reflection and 
transmission coefficients for graphene. These coefficients can 
be found by matching the magnetic H and electric E fields 
through the boundary conditions: 



e z x (E_ - E+) = 0, 

47T 47T (7) 

e z x (H_ - H+) = — j = — -ae z x [e z x E+], 



where E_ (H_) and E + (H + ) stay for the electric (magnetic) 
fields in the regions of negative and positive z, respectively, 
and e z is the unitary vector along the +z direction. As a result 
of the matching we have 



TE 
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nTE Qz 
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TM 
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aq z + 1 ' 
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(8) 



with a = 27r<j/c, being the dimensionless 2D conductivity. 
Explicitly, we have for G = Gq M + Gq E 



8^2 
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QzQ 2 



ic\r+iq z \ Z - Z '\ 



( $ 



-q x q y 
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G™(r) 
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/ QxQz q x q y q z Tq x q 2 



q x q y q z q y q z Tq y q 

9 ri 2 ^4 , 



\Tq x q Tq y q q /q z 



Analogously, G R = G R + Gjf reads 



f r>TE iqr-iq z (z'+z) 

8* 2 J q z q 2 q 

( q 2 y -q x q y o 
2 o 



V 



-q x q y 
o 



q x 

o o ; 
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dq 



pTM ic l r-iq z ( Z / -\-z) 
2 Q 



(10) 



8tt 2 J q 

( q 2 x q z q x q y q z -q x q 2 



q x q y q z q y q z 



\ q x q 
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-q y q 
-q A /q z 



and finally, for the transmission part Gt = G T + G T we 
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have 



G T T E (v) = 



ik u . 



I 



rfq 



8tt 2 J q z q 2 q 



rpTE giqr+iq z {z'+z) 



( q% —QxQy o\ 



V 



G? M (r) 



ikt. 



IxQy ql 
; 

dtlrpTM iqr+iq z (z'+z) 

Sn 2 J q 2 q 

( q 2 x q z q x q y q z -q x q 2 

q x q y q z q 2 y q z ~q y q 2 

q A /q z 



(ID 



q 2 y q z 

-QyQ 2 



\-q x q 

Eqs. (|9])-([TT]) present the angular representation of DGF in 
Cartesian coordinates. By transforming a cylindrical coordi- 
nate system, these expression can be greatly simplified. 




Fig. 2. Transformation of the coordinates. 



B. DGF in cylindrical coordinates 

As was previously shown, the Purcell factor (the total decay 
rate normalized to the free space decay rate) of the point 
emitter placed directly over the graphene monolayer diverges 
due to losses through the evanescent waves with large q- 
components (6), (7), fT0| . This means that the dipole placed 
directly on the graphene surface is quenched. However, the 
parameter that accounts for the efficiency of the coupling to 
GSP, (J}-f actor, defined as the ratio of the emitter's decay 
rate through GSP to its total decay rate) has an optimum 
value in the region of very small distances from the dipole 
to graphene: \z f opt \/X ~ 10 -2 (see [ 10 1). Then, taking into 
account that the problems related to the high values of /3-factor 
are relevant, the range of small distances presents a special 
interest. Another important point is that the field patterns at 
the distances r > \z'\ do not differ essentially upon the field 
patterns created by a dipole lying directly on the monolayer. 



Therefore, for the above two reasons, in this paper we will 
consider the dipole placed directly onto the graphene sheet. 

We would like to notice, that some physical systems can be 
reduced to a problem of a dipole lying directly on the graphene 
sheet. For instance, a subwavelength aperture in graphene 
sheet can be represented by an effective dipole placed directly 
onto the sheet, (for comparison with the case of metal films 
see (2TJ, (22)). 

Additionally, since the analytical treatment of both reflection 
(z < 0) and transmission (z > 0) parts of the DGF is similar, 
we will derive the expressions for the transmission part, Gt- 

So, supposing that z' = 0~ in the previous expressions, the 
Green's dyadic G(r, z) = Gt{y, z, z' = 0~) can be simplified. 
The symmetry of the problem makes it convenient to work in 
cylindrical coordinates (r, 0, z), see Fig. [2] (a): 



r cos ( 



y 



r sin ( 



z. 



(12) 



In this system of coordinates the Green's dyadic can be 
obtained from the one in cartesian coordinates through 



Q°yl rp—l Qcart rp 



where 



/C0S( 

sin c 

V 



— sm ( 
cos <j) 




(13) 



(14) 



Having performed this transformation, the DGF in cylindrical 
coordinates is expressed in the form of Sommerfeld integrals 

as G — Gr> -\- G s \ 



G p (r,z) 
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(15) 
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where the subscripts "s" and "p" correspond to TE and 
TM polarizations respectively. In this expressions, J±(qr) = 
Jo(qr)±J2(qr) and J n (qr) are Bessel functions of nth order. 
Equations ( [T5] ) can be treated numerically. We show some 
recipes for the integration in the complex plane q in Appendix 
A. Let us now proceed with the asymptotic expansion of the 
DGF. 

IV. Asymptotic expansion of the DGF 

In this section we will derive explicit asymptotic expressions 
for the elements of DGF following the steepest-decent method, 
modified in order to take into account the presence of both 
poles and branch points close to the integration path (13), 

CD- 
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First, using the identity 

2J n = HV>+HV>, 



(16) 



i[x-f (n+i)] I l + i 



An 2 



8x 



+ 0{x-i). 
(18) 



For convenience, let us normalize the DGF as follows 

AT 



G 



8tt 



— 9- 



(19) 



Then using ([T6j)-(|T7]), we find from ( p3] ) the expressions for 



9 = 9s+g T 
9r(r, z) 



dq 



3 iqr+iq z z 



(20) 



with r = s,p and the superscripts (1, 2) of A indicate distinct 
dependencies upon r. The denominators in the integral are 



f s (q) =a + q z , f p (q) = aq z + 1, 
and nominators 



(21) 



q z rr — q (rz + zr) H zz 



1 



i(2) 



ii 2 ) = 



lq z rr — Sq 
-rr H — 00 



3g (rz + if) z£ 



(22) 



where r , 0, z represent the unit vectors. 

As we see, the integrands contain branch points correspond- 
ing to q z , and branch cuts corresponding to lm(q z ) = 0. In 
order to remove these problematic peculiarities, we perform 
the following standard change of integration variable and 
coordinates (see Fig. [2](b)): 



where H n are Hankel functions, we can extend the limit of 
integration to the whole real g-axis in ( [T5] ) 

pOC -i /'CO 9" 

/ dqJ n (qr)F(q) = - dq (qr)F(q), (17) | 

where according to Eqs. ([15) function F(q) is odd/even for 
even/odd values of n. In (fTT] ) we have used the symmetry of 
the Hankel functions H^^J-x) = -e in7T (x). 

Second, we use the asymptotic form of the Hankel functions 
for large arguments (qr ^> 1). Notice that the region that 
provides the major contribution to the integral corresponds to 
q > 1. Then the formal condition of the asymptotic expansion 
validity for the lower value of the contributing q reads as 
r ^> 1. However, as we will show below, the true region of 
distances where the asymptotic approximation is valid is much 
less restricted. Retaining the first two terms in Hn\qr), the ^ 
expansion reads w 
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Fig. 3. Integration contours in the complex plane of the integration variables 
(p [panel (a)] and w [panel (b)] for different values of angle 0. In panels (a) 
and (b) the initial integration contour corresponds to C = C\ + C2 + C3 



and C = C' 



o 2 



C 3 



respectively. The pole position is shown by a 



circular symbol. In the complex plane of (p the initial path is the same and 
the steepest-decent path depends upon 0, while in the complex plane w the 
steepest-decent path is the same, lm(w) = 0, and the initial integration path 
changes with 0. The position of the pole in the complex plane w is dependent 
upon 6. 



The integrals ([20]) then transforms as 

00 

j (]<!<' 



A M (a) 



fr(q) 



(24) 



dcpe 



ip cos(ip — 9) 



cos ip 



i^ n) (sm(^) 



ic /r(sin^) 

where the integration contour C pases through the complex 
plane ip, see Fig. [3] (a) and corresponds to the real axis 
in the complex g-plane. At this stage we can perform a 
steepest-decent integration. Then the integration path must be 
transformed to cos [Re (9?) — 0] cosh [Im ((p)] = 1, see Fig. [3] (a), 
with the saddle point cp = 6. However, the steepest-decent 
integration is much easier in another complex variable plane. 
This variable w is given as follows 



w 



V2e % 4 sin 



(25) 



sm ip, 



cos ip, 



ps'mO, z = p cos 6. (23) so that i cos((p — 6) = i — w 2 . With this change of variable g T 
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becomes 



= e ip [ 
Jc 



dwe 



-pw 



r 



dtp 
dw 



(26) 



fr[q(w)] ' 

where the integration contour C is shown in see Fig. [3] (b). 
The steepest-decent integration path is now simply given by 
\m.(w) = 0, and the saddle point is located in the origin, 
w = 0. Notice that the branch points q = ±1 (corresponding to 
q z = 0) are given by tp = ±tt/2 in the complex plane tp, while 
in the w-plane they are located at w = \/2e 1 ^ sin (± \ — |). 
For 6 = ±7r/2 the branch point and saddle point in u>-plane 
coincide. 

An important point here is that the elements of the integrant 
dyadic in g T are singular due to the presence of the poles 
Up(o) = an d f s (o) — 0] m me denominators. These poles 
are located at q = q p and q = q s respectively. They correspond 
to TM-surface wave (GSP) and TE surface wave (2) 



1 



(27) 



The positions of the poles depend essentially upon the value 
of the normalized conductivity a. TM waves correspond to 
Im(a) > 0, while TE ones correspond to Im(a) < 0, so 
that TM and TE surface waves cannot exist at the same 
frequency see [2]. High values of \a\ correspond to large q s 
while low values of \a\ correspond to large q p . In principle, 
taking into account a wide range of metamaterials that are 
available at present, a wide range of a is also accessible. 
A mathematical treatment of the problem corresponding to 
a three-dimensional (3D) layer of a very thin thickness h <C A 
with the dielectric permittivity e^d can be reduced to the case 
of a 2D sheet with an effective 2D normalized conductivity 
a e ff. The relation between 2D effective conductivity and 3D 
permittivity is established from the comparison of the Fresnel 
coefficients and reads as a e ff = Trheso/i^- 

Returning to the case of graphene, due to small values of 
\a\ in graphene, TE surface waves are very weakly bounded 
(\q s \ ~ 1) and therefore they virtually do not couple to the 
point emitter. This means that, in practice, the contribution 
from TE pole can be neglected in graphene |8|. Nevertheless, 
taking into account a wide range of possible 2D sheets, in our 
analysis we retain the contribution from both poles. 

In order to proceed with the series expansion, the singular 
terms in the integrands must be separated. The separation for 
the dyadics ^^(w) into a pole and a smooth part &^( w ) 
is as follows: 



Q 



(n) 



w — W T 



(28) 



w — W T 



" (n) 

where Q K T ' are dyadics with the elements corresponding to 
the residues of ^^(w), which after some algebra can be 



computed from Eq. (|26[) as: 
1 



4 n) fe), Qi n) = -ii n) fe) 



(29) 



a*q p ' q s 

Let us explicitly write out the expressions for the residue 
dyadics 



0(1) _ X^L 



Q 



(1) 



Qp 



(2) 



(2) 




(30) 



Q 



Now we can deform the integration contour C in the w 
complex plane, into the real w axis. Then the singular terms 



Qt 



in Eq. ([28]) for &r ,2 \w) that enter to the integral ([28]) 
can be integrated analytically. The result for g T (r,z) can be 
presented in the form of a sum 



g T (r, z) 



in e 



,p(i-wl) 



erfc(—iw ry /p) ( Q 



)(2) 



+ <?0t0, Z), 



(31) 

where the term with the complementary error function 
erfc(x) = (2/y/ir) e~ l dt is due the singularity. This 
function includes the contribution from the pole that is auto- 
matically taken into account when the pole is crossed by the 
transformation of the integration contour. The second term in 
Eq. ( [3T] ), go T (r,z) presents a nonsingular contribution 



).«»/ 



dw e~ pw 



(32) 



The integral appearing in Eq. p2| ) is of the Gauss type, so the 
functions 4>^:(w) can be expanded in Tailor series close to 
w = and integration of every term is easily performed with 
the following result 

i r(*±g) 

(33) 



9or(r,z 
d 



dw n 



w=0 



where T is Gamma function. Retaining in this expression the 
terms up to order r -3 / 2 (which is enough for the most of the 
cases) yields 



gor(r, z) 



(0) 



Ap dw 2 



(34) 
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Recalling that the saddle point w = corresponds to ip = 9 
and therefore to q = sin 6, we can explicitly write out the 
dyadics 6^(0) 



*W(0) = 



V2e 



,ijg(sing) 
/r(sine) 



cost 



(35) 



where we have taken into account the identity ^ = 
V2e~ 1 ^ I cos(^^). The explicit expression for the third term 
in ( [34] ) is more cumbersome for arbitrary 6, therefore we 
give its formal expression, involving derivatives of previously 
defined functions 



dw 2 
dw 2 



2Q« 



\w=0 



=0 2v^e" 



d 2 ^\ 
dw 2 

4 d^ 2 



|w=Ch 



cos (sin ip) 

cos(^) /r(sin^) 
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The expressions pT) , with g 0r (r,z) given by ( [34] ), present 
the analytical approximation of the DGF. Let us now analyze 
different terms in this expression and the applicability of the 
approximation. 



When the graphene sheet disappears, a — ^ 0, we recover 
the spherical wave term (~ 1/p) of DGF corresponding to 
a dipole in FS. Notice that at 9 = tt/2 (for the fields along 
z = 0), the elements rr, rz, zr in Gp° vanish independently 
upon whether the graphene sheet is present or not. In contrast, 
the element (jxj) in Gf° at 9 = ±tt/2 is very sensitive to the 
presence of graphene. It takes non-zero values for free space, 
a = and vanishes for a ^ 0. This property of the (jxj) 
element is similar to the diffraction shadow effect in metals 
due to the presence of surface modes and formation of the 
Norton waves [19 ]. In case of metals, however, the diffraction 
shadow appears for the TM part of the dyadic, while in thin 
films this takes place for the TE part. 

The fact that some elements of the dyadic G RO vanish 
indicates that other terms in the general solution for G 
must be considered in order to provide the correct far-field 
representation of the DGF and electric fields. Let us consider 
in details the case of 9 = tt/2. 



A. Asymptotic expansion of DGF in graphene plane, z = 
0- (9 = tt/2) 



V. Analysis of the analytical solution 

In its general form, the analytical solution presents a non- 
trivial combination of the contribution from the pole and 
saddle point. The "interaction" between these contributions 
depends both upon the distance between the saddle point and 
the pole in the complex plane and upon the physical distance r 
responsible for the oscillations of the integrand. The parameter 
that measures the interaction between the saddle point and the 
pole is called by Sommerfeld "numerical distance" d T and 
its square presents the argument of the complementary error 



function, d 
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In the g-plane, the saddle point corresponds to the condition 
of the extremum of the exponential phase in pQ| , i.e. to 
(Q/Qz)min — r/z, or q m in = sin#. This can be considered 
as an equation for the rays in "Ray Optics" (RO). If the 
contribution of the pole is neglected, then the same result can 
be derived following the standard stationary phase evaluation. 
In the far field, the leading RO contribution corresponds to the 



first term in ([34]) with <&^(0) replaced by $V ±; (0). It reads 



9?°(r,z) 



/2tt i^ n) (sin^ 



f ( m cos ^ (3?) 
p / r (sin0) 

Returning to the DGF via Eq. ( fT9| ), we have explicitly [using 
Eqs. ( |2T] ), ( [22] ) and the relation r = p sin 0] 

( cosO — sin#\ 

k UJ e' lp cosO 



47rp(l + a cos 6) 



k 0J e lp cos 9 
47rp(a + cos 6 



V" 
/o 








sin 2 



(38) 



Here we present simplified expressions for g T (r, z) at z = 0. 
Recall that the general expression is given by the Eq. ( [3T] ), with 
go T given by Eq. ([34]). The functions 4>r n ^(0) at 9 = tt/2 read 



$W(0) = $f(0) = 0, 

4>«(0) = yfar^zz, $< i ? ) (0) = -^e- i Tzz. ^ 



Introducing for a brevity of notations the following dyadic 



dw 2 1 



=0 = M T for 



(40) 



The derivatives <|36]> simplify as 
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Now we can explicitly write the full expression for the TE 
dyadic 
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(43) 



Recall that in Eqs. ([42|)-([43]) the residue dyadics are given 
by Eq. ( [30| ), and M r have the form of Eq. ( |4T] ). The locations 
of the poles in the complex w-plane at 6 = tt/2 read 



w s 



e ^V^-l. 



(44) 



We have checked that the Eqs. ( |42| ), ( [43] ) transform to the 
DGF of FS in the limit a — )> 0. To perform this limit, one 
should carefully expand all the coefficients taking into account 



that for small a the poles become w s ~ ^ e * 4 ' 



y/a' 



In particular, for TM part of the DGF the expansion for large 
arguments of the complementary error function must be taken. 

B. Numerical check on the validity of the analytical approxi- 
mation 

Let us present some illustrative examples that demonstrate 
the validity of the analytical approximation. For this we 
directly compare numerical and analytical calculations for two 
components of DGF in the most unfavorable situation, i.e. for 
z = (6 = 7r/2), when the "RO" contribution disappears. 
We perform the precise (converged) numerical calculations 
according to Appendix A. The analytical approximation is 
given by the expressions ( |42| ), ( [43] ). 

The conductivity of graphene is a function of frequency, v = 
oj/(2tt), chemical potential p, temperature T and scattering 
time r, see Appendix B. For the illustration, we consider the 
room temperature T = 300K and p = 0.2eV (48THz), which 
are typical experimental values. 

The comparison is shown in Fig. [4] The range of the 
distances corresponds to the subwavelength region (outside 
of this region the the numerical and analytical curves are 
virtually undistinguishable). In order to characterize the dif- 
ference between numerical, G n , and analytical, G a , results, 
in the insets we have represented the relative error A^/ = 



in THz frequency range. According 
to the insets to Fig. [4j the error is a non-monotonous function 
of frequency. However, it has a decaying tendency with 
frequency increase. To understand such behavior, let us notice 
that for higher frequencies \a\ decreases so that both the real 
and imaginary parts of q p increase. In the lower limit (y = 1 
THz) a ~ 0.11 + 0.69i so that q p ~ 1.7 + 0.19i, while in 
the upper limit (y = 10 THz) a ~ 0.0016 + 0.07z so that 
q p ~ 14.34 + 0.34z. Thus, the propagation length of the GSP, 
Lgsp = A/[27rIm(#p)], decreases due to increase of Im(q p ). 
This leads to a strong spacial decay of the terms in the solution 
(at the deep sub- wavelength distance), related to the GSP field 
components. Since the analytical solution recovers the FS DGF 




Fig. 4. Comparison between the numeric and analytic calculations of zz and 
zr DGF elements at z = 0. The main figures show both real and imaginary 
parts of G zz and G zr as a function of distance. The values of these elements 
have been normalized to the maximal values of their modules in the shown 
range of distances. The parameters for graphene in both panels are: T = 
300K, /i = 0.2 eV, r = lps, v = 10 THz. The insets show the dependencies 
of the relative error upon the frequency for different distances. In the insets 
\i and r are the same as in main figures. 



(up to 1/r 2 ), the coincidence between the analytical and exact 
solutions improves for higher frequencies. 

VI. Long distance limit. Surface modes and 

ALGEBRAICALLY-DECAYING COMPONENTS 

In the region of parameters, where the argument of the 
complementary error function is a large number \w T \yfp ^> 1 
(the numerical distance is large, ^> 1), this function can 
be substituted by a few terms from its asymptotic expansion 



erfc(— iw T y/p) = 20 [Im(it; r )] 
e w *P g (-l) n (2n)\ 



(45) 



'^V^P (-0 2n+1 n\(2w T ) 2n p n ' 

9_O) = 0, for £>0, 
6_(x) = l, for £<0. 



where 



(46) 



The first term (which is independent upon p), appears when 
the transformation of the initial integration path to the steepest 
descent one results in crossing the pole. Retaining the terms 
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exact up to r 3 / 2 in Eq. ([45}, we arrive at 

erfc(— iw r y/~p) = 20_ [Im(iu r )] -f 
ie w - p ( 1 



1 



■o[KVp)" 



(47) 



Let us concentrate on the case of the in-plane fields, 6 — it/ 2. 
Taking into account that e w T r e l( i^ r = e ir ? the expressions ( |42} , 
(|43} simplify to 



G s (r,0) 



fc w e*4 
8tt 



G p (r,0) = 



n, e" 



4ry / r 



M, 



1 

(48) 




4ry/r 



o 4 



(49) 



where the dyadic II r describes the surface mode and only 
contributes when the pole is located in the physically proper 
Riemann sheet 



n r = 2iri • 9_ [lm(w T )] QV 



)(!) 



)(2) 



(50) 



We would like to notice that in the approximate expression 
for g s (r, 0) given by Eq. ( |48] > we cannot recover the limit 
a — )> anymore. 

Expressions ( [49] ), ( [50] ) present a sum of the surface mode 
term proportional to II r (GSP in case of p-polarization) and 
algebraically-decaying terms. The term describing the surface 
mode can be directly recovered from the angular representation 
Eq. ( [20] ), considering only the residue of the pole. The alge- 
braic components, in their turn can be derived from the same 
integral p0| , considering the contribution from the branch- 
point q z — (see Ref. (81). Asymptotic expressions ( |49} , ( [50] ) 
present thus independent contributions from the pole and the 
branch cut. 

The main physical reason of the validity of this approx- 
imation is that for sufficiently long distances only sharp 
peculiarities on the density of electromagnetic states (DES) 
contribute. DES is reflected by the integrand in Eq. ([20}. For 
large r, the smooth region of (DES) is progressively canceled 
out in the integral, which is eventually dominated by the strong 
(and rapid) contribution from the pole. The contribution of this 
pole gives the field of the surface mode. Due to losses, the 
density of states associated with the pole has a finite width, 
which causes the exponential decrease of the GSP amplitude 
with distance (characterized by the surface mode propagation 
length). Then, the contribution to the integral from either 
kink or square-root singularity (ex l/q z ) located at q z = 
dominates. This takes place since this kind of features cannot 
be characterized by a typical width in q- space (these features 
are infinitely sharp in q- space), and they are not as strongly 
suppressed as the pole contribution when integrated with an 
oscillatory function. The contribution of the kink/square-root 
singularity yields the algebraic decay of the DGF components 
with respect to the distance. 



A detailed physical description of the algebraically- 
decaying components of the fields from a point source can be 
found in Ref. (8). Let us recall here the physical meaning of 
all the algebraically-decaying terms in Eqs. ([48}, ( [50} . These 
terms appear in the dyadics M T , and the element zz of G p , 
contains an additional contribution ~ 1/r. The algebraically 
decaying components are composed of both FS terms and 
Norton waves (NW) [ 20 1. The FS terms do not depend upon a 
and result from the contribution of the branch-cut singularity 
(l/q z ) that yields the dependency ~ 1/r and a kink that 
yields ~ 1/r 2 . Notice that while zz component (of the TM 
part) contains both the term decaying as ~ 1/r and 1/r 2 , 
the FS part of the element rr has only 1/r 2 decay. All the 
rest of the 1/r 2 terms that depend upon a correspond to the 
NW (compare with the case of metals, Ref. [19], (2TJ). We 
would like to notice that as follows from Eqs. ( [48} (where a 
has been supposed to have a nonzero value) the element (jxj) 
contains only the NW. However, if we carefully perform the 
limit a — » in the initial equation ([42}, recovering the FS 
DGF, the element (jxj) will contain a ~ 1/r term. This can 
be explained by the fact that in the DGF given by its angular 
representation ( [20} has a square-root singularity for a = 0. 

Since the main message of this paper is the analytical 
treatment of DGF for graphene, let us illustrate the validity 
of the asymptotic expressions ([49}, ( [50} by comparing two 
elements (zr and zz) of DGF with the numeric solution. 

First, the competition between algebraically-decaying and 
GSP terms is shown in Fig. [5] (a,b). At the beginning of the 
shown spacial window, both elements of DGF are dominated 
by the GSP terms. Then, in the region of R ~ (7 — 9) A 
for G zr and R ~ (4 — 5) A there is a crossover, where the 
exponentially decaying GSP is overcome by the algebraically 
decaying field. The algebraic decay for G zr corresponds only 
to NW (~ l/R 2 ), since the FS contribution is zero for 
this element. In contrast, the asymptotic behaviour for G zz 
corresponds to FS with dominating ~ 1/R term. The NW term 
is also present in the element zz, but its contribution is much 
weaker than that of the FS component. In the region of the 
cross-over the field possesses a peculiar two-scaled oscillation 
behavior (corresponding to the wavelength of the GSP and 
vacuum wavelength). Notice that the amplitude of the field 
at the crossover is extremely small, so for this instance the 
analysis has mainly an academic value. 

Second, in Fig. [5] (c,d) a direct comparison of numeric and 
asymptotic results is performed in the interference region. As 
one can see, the asymptotic approximation perfectly captures 
all the details of the exact result. We have also checked the 
validity of our asymptotical expressions of all other DGF 
elements. 

VII. Conclusion 

We have performed an analytical treatment of Dyadic 
Green's Function for 2D sheet. In particular, we have tested 
the analytical expressions on the case of graphene. We have 
retained all the necessary terms that provide high precision 
(~ 1%) down to distances of 1/10 wavelengths and reasonable 
precision (~ 10%) down to 1/100 wavelength. 
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Fig. 5. (a,b) The absolute values of G zr and G zz as a function of the distance from the point source. The modulus of GSP, NW in (a) and GSP, NW in 
(b) terms are also rendered in the same panels. The insets to (a,b) present zooms of the main panels in the region of a strong interference between GSP and 
algebraically-decaying components (FS and NW). (c,d) The real part of G zr and G zz as a function of the distance from the point source in the region of 
interference between GSP and NW or GSP and FS. The analytical expressions based upon Eq. p9) are compared with the exact calculation. In both (a,b) 
and (c,d) all represented values are normalized to the maximal value of \G zr \ (\G ZZ \) in the shown intervals. The parameters of graphene are the same as 
in Fig. g 



For the limit of long distances (in units of plasmon wave- 
lengths) we have presented simplified expressions with sepa- 
rated contribution from the pole (plasmonic field) and from the 
branch point (algebraically decaying field). These expressions 
are relevant for future studies of the electromagnetic properties 
of sub wavelengths objects placed on a graphene sheet. 

Appendix A 

Numerical computations of Sommerfeld integrals 
Let us consider the integral of the following form 



J(r) 



dqF(q, r), F(q, r) = F(q)J n (qr), (51) 



where J n is the nth-order Bessel function and the function 
F(q) remains finite for Im(g) — >> oo. We suppose that the 
function F(q) has a pole at q = q p and is dependent upon 
q z = y/l — q 2 so that it has branch cuts Im(q z ) = 0. The 
pole and the branch cut are not the only difficulties of the 
integral. In case of the integration along the real axis of the 
complex g-plane, the Bessel function has a strong oscillatory 
behavior for q ^> 1 and the integration is very delicate. When 
r increases, the convergence of the integral becomes worse. 
In order to stay away from the singularity and remain at the 
same Riemann sheet, the integration path can be deformed 



according to Cauchy theorem (supposing that we do not cross 
the pole) 



J(r) = J dqF(q)J n (qr) + ±J dqF(q)H^(qr) 

+ \ f dqF{q)H^{qr). 
J c 



(52) 



Here the contour A passes below the real axis rounding the 
pole and the branch cut and then returns towards the real axis 
at the point q = S with S > Re(q p ), moving into lm(q) — >> oo 
for Hn^ term, and to lm(q) — » — oo for (see Fig. [5] the 
paths marked by "2"). The contours "B" and "C" are restricted 
by the limiting values lm(q) — ±A. 

When the pole is far away from the origin, \q v \ ^> 1, or/and 
the distance parameter is large, r ^> 1, it is convenient to 
bend the contours before the pole, i.e. choose S < Re(q p ) 
(see Fig. [6j the paths marked by "1"). In this case the pole 
in the second integral of Eq. ( |52| ) must be taken into account. 
The contribution of the pole adds the residue term into the 
expression ([52|): 



J(r) 



J A dq... + \j B dq... + \j cd q... 

+iri-Res(F,q p )H ( n 1) (q p r). 



(53) 
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where the intraband and interband contributions are: 



-2 



Fig. 6. The contours for the integrand corresponding to the DGF elements 
G rz = G zr (the real part of the TM-term) for z = 0, r = 0.5. The integrand 
is normalized to the maximal value of its module. In the region lm(q) > 
(lm(q) < 0) the Bessel function J\ is replaced by the decaying Hankel 
function (H^). As a result of this replacement, a discontinuity along 
lm(q) = appears. The pole position q p = 2 + 0.5i corresponds to a 
"toy value" of the normalized conductivity a ~ 0.168 + 0.516i, chosen 
for better visualization. For the same reasons, in order to better il lustrate the 
branch cut Im(q z ) = 0, the casuality has been exaggerated: \ — q 2 — > 

a/(T+ iO) 2 — q 2 — > y/(l + O.H) 2 — q 2 . The parameters of the contours: 
8i = 1.5, 5 2 = 2.75; wi = 1, w 2 = 1.5; Ai = A 2 = 2. 



Each path "A", "B" and "C" can be parameterized qi = qi(t) 
(i = A, £?, C) so that the integration is reduced to the domain 

[0,1]: 

1 

dt 



i 

J2 J dtF[ qi {t)} [ 



with 



QA(t) = 6 • t — iw sin (nt) . 
q B (t) = 5 + itA, 
qc(t) = 6-itA, 



(55) 



where the parameters 5, w and A are chosen so that the best 
convergency of the integrals is provided. In this paper, we have 
performed the integration over t following Simpson's rule. 

Appendix B 
Graphene's conductivity model 

The conductivity of graphene computed within the random 
phase approximation |T6|-p8) can be written through the 
chemical potential /i, the temperature T, and the scattering 
energy £ s as follows 



2ieH 



O 'intra 



&inte 



2tt ln (n 



•27) 
1 



ln 



arctan 



2 cosh 

n 



■2) 



2t 



(57) 



2) 2 + (2t) 2 



In this expressions Q = huj/jjL, 7 = £ 8 /fi and t = T/fi, with 
T expressed in units of energy. The scattering energy is related 
to the relaxation time r as r = £ s /%. 
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